The Imprint of The Extragalactic Background Light in 
the Gamma-Ray Spectra of Blazars 
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The light emitted by stars and accreting compact objects through the history 
of the Universe is encoded in the intensity of the extragalactic background light 
(EBL). Knowledge of the EBL is important to understand the nature of star 
formation and galaxy evolution, but direct measurements of the EBL are lim- 
ited by Galactic and other foreground emissions. Here we report an absorption 
feature seen in the combined spectra of a sample of gamma-ray blazars out to 
a redshift of z~1.6. This feature is caused by attenuation of gamma rays by 
the EBL at optical to UV frequencies, and allowed us to measure the EBL flux 
density in this frequency band. 

The bulk of the intergalactic gas in the Universe must have been reionized between the 
epoch of cosmic recombination, when the Universe was only 300,000 years old (z~l 100), and 
1 billion years later (z~6) as indicated observationally by the spectra of distant quasi-stellar 
objects (1). However, the sources, modes and nature of this cosmic reionization are largely 
unknown because most of this redshift range has yet to be explored. Photoionization by UV 
radiation, produced by the first stars and galaxies of the Universe, represents the primary suspect 
for the ionizing process (2, 3). Direct detection of the UV radiation fields is thus of fundamental 
importance, but at present extremely difficult (3). 

An indirect but powerful means of probing the diffuse radiation fields is through 7-7 ab- 
sorption of high-energy gamma rays (4-6). In this process, a gamma-ray photon of energy E 1 
and an EBL photon of energy E EBL annihilate and create an electron-positron pair. This pro- 
cess occurs for head-on collisions when (e.g.) E 1 x E EBL > 2(m e c 2 ) 2 , where m e c 2 is the rest 
mass energy of the electron. This introduces an attenuation in the spectra of gamma-ray sources 
above a critical gamma-ray energy of E crit (z) ~ 170(1 + z)~ 2 - 38 GeV (7, 8). 

The detection of the gamma-ray horizon (i.e. the point beyond which the emission of 
gamma-ray sources is strongly attenuated) is one of the primary scientific drivers of the Fermi 
Gamma-ray Space Telescope (9-11). Several attempts have been made in the past but none 
detected the long-sought EBL attenuation (12-14). So far, limits on the EBL density have been 
inferred from the absence of absorption features in the spectra of individual blazars (13, 15), 
distant galaxies with bright gamma-ray emission powered by matter accreting onto central, 
massive black holes. While this feature is indeed difficult to constrain for a single source, we 
show that it is detected collectively in the gamma-ray spectra of a sample of blazars as a cut-off 
that changes amplitude and energy with redshift. We searched for an attenuation of the spectra 
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of blazars in the 1-500 GeV band using the first 46 months of observations of the Large Area 
Telescope (LAT) on board the Fermi satellite. At these energies gamma rays are absorbed by 
EBL photons in the optical to UV range. Thanks to the large energy and redshift coverage, 
Fermi-LAT measures the intrinsic (i.e. unabsorbed) spectrum up to ~ 100 GeV for any blazar at 
z<0.2, and up to ~15 GeV for any redshift. 

The LAT has detected > 1000 blazars to date (16). We restricted our search to a subset of 150 
blazars of the BL Lacertae (BL Lac) type that are significantly detected above 3 GeV, because 
of the expected lack of intrinsic absorption (17). The sample covers a redshift range 0.03- 
1.6 (18, 19). The critical energy is therefore always >25 GeV, which means that the spectrum 
measured below this energy is unabsorbed and a true representation of the intrinsic spectrum of 
the source. We thus determined the intrinsic source spectrum relying on data between 1 GeV 
and the critical energy E crit and extrapolated it to higher energies. By combining all the spectra 
we were able to determine, the average deviation, above the critical energy, of the measured 
spectra from the intrinsic ones, which ultimately provides a measurement of the optical depth 

T 77 . 

The analysis was performed using the Fermi Science Tools (20). We determined the spectral 
parameters of each blazar by maximizing the likelihood of a given source model. The model 
comprised the Galactic and isotropic diffuse components and all sources in the second Fermi 
LAT catalog (21) within a region of interest (ROI) of 15° radius. We modeled the spectra of the 
sources in our sample as parabolic in the logarithmic space of energy and flux (see Eq. 2 in (21) 
for a definition). Their spectra were modified by a term e~ T ^^ E,z ' that describes the absorption 
of gamma-ray photons on the EBL. In the above we defined t 7J (E,z) = b • r™ ode '(_E, z), 
where the r™° del (E,z) is the optical depth predicted by EBL models (7, 22-25) and b is a 
scaling variable, left free in the likelihood maximization. In particular, this allowed us to assess 
the likelihood of two important scenarios: i) there is no EBL attenuation (6=0), ii) the model 
prediction is correct (6=1). 

We combined the data from all the ROIs in a global fit that determined the common param- 
eter b for a given EBL model (see Table SI). All those models with a minimal EBL density 
based on (or compatible with) resolved galaxy counts (2, 7, 24-27) were found to be acceptable 
descriptions of the Fermi data (i.e are consistent with b=\ within ^25 %, see also Figured]) 
yielding a significance of the absorption feature of up to ~6 a. Models that predict a larger 
intensity of the EBL particularly in the UV (22, 23) would produce a stronger-than-observed 
attenuation feature and are therefore incompatible with the Fermi observations. Our measure- 
ment points to a minimal level of the optical-UV EBL up to redshift z^l.6 which combined 
with the upper limits (15, 28, 29) derived at lower redshift (using observations of blazars at TeV 
energies) on the near-infrared EBL highlights the conclusion that most of the EBL intensity can 
be explained by the measured galaxy emission. 

Our measurement relies on the accuracy of the extrapolation of the intrinsic spectra of the 
sources above the critical energy (30). This in turn depends on a precise description of the 
gamma-ray spectra by our source parametrization. To verify that this is the case and to exclude 
the possibility that the detected absorption feature is intrinsic to the gamma-ray sources (17), we 
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performed the analysis in 3 independent redshift intervals (z<0.2, 0.2<z<0.5, and 0.5<z<1.6). 
The deviations from the intrinsic spectra in the three redshift intervals are displayed in Figured 
In the local Universe (z<0.2), EBL absorption is negligible in most of the Fermi-LAT energy 
band (E crit > 120 GeV). The lowest redshift interval therefore reveals directly the intrinsic spec- 
tra of the sources and shows that our spectral parametrization is accurate (18). The absorption 
feature is clearly visible above the critical energy in the higher redshift bins. Its amplitude 
and modulation in energy evolve with redshift as expected for EBL absorption. In principle, 
the observed attenuation could be due to a spectral cut-off that is intrinsic to the gamma-ray 
sources. The absence of a cut-off in the spectra of sources with z<0.2 would require that the 
properties of BL Lacs change with redshift or luminosity. It remains an issue of debate whether 
such evolution exists (31-34). However, in case it were present, the intrinsic cut-off would be 
expected to evolve differently with redshift than we observe. To illustrate this effect, we fitted 
the blazar sample assuming that all the sources have an exponential cut-off at an energy E . 
From source to source the observed cut-off energy changes because of the source redshift and 
because we assumed that blazars as a population are distributed in a sequence such as that pro- 
posed in (31-34). E was fitted to the data globally like b above. As apparent from Figured it 
appears difficult to reconcile the observed feature with an intrinsic characteristic of the blazars' 
spectra. We therefore associate the spectral feature to the EBL absorption. 

At energies < 100 GeV, gamma rays observed at Earth and coming from redshift >1 inter- 
act mostly with UV photons of >5 electron volts. An UV background in excess of the light 
emitted by resolved galaxies can be produced locally by AGN or at higher redshift (z^7-15) by 
low-metallicity massive stars (35). By comparing the results from the best-fit EBL models, we 
measured the UV component of the EBL to have an intensity of 3(±l)nW m~ 2 sr _1 at z«l. 
A contribution to the UV background from AGN as large as the one predicted by (36) (i.e., 
ps 10 nW m~ 2 sr -1 ) and used in the EBL model of (22) is thus excluded by our analysis at high 
confidence. However, the recent prediction (37) of the UV background from AGN (ps 2 nW 
m~ 2 sr -1 ) is in agreement with the Fermi measurement. Direct measurements of the extra- 
galactic UV background are hampered by the strong dust-scattered Galactic radiation (38). The 
agreement between the intensity of the UV background as measured with Fermi and that due to 
galaxies individually resolved by the Hubble Space Telescope (39) (3±1 nW m~ 2 sr -1 versus 
2.9-3.9 nW m~ 2 sr -1 , respectively) shows that the room for any residual diffuse UV emission is 
small. This conclusion is reinforced by the good agreement of the Fermi measurement and the 
estimate of the average UV background, at z>1.7, of 2.2-4.0 nW m~ 2 sr -1 using the proximity 
effect in quasar spectra (40). 

Zero-metallicity population-Ill stars or low-metallicity population-II stars are thought to be 
the first stars to form in the Universe and formally marked the end of the dark ages when, with 
their UV light, these objects started ionizing the intergalactic medium (41). These stars, whose 
mass might have exceeded one hundred times the mass of our Sun, are also believed to be 
responsible for creating the first metals and dispersing them in the intergalactic medium (42- 
44). A very large contribution of population-Ill stars to the near-infrared EBL had already been 
excluded by (15). Our measurement constrains, according to (45, 46), the redshift of maximum 



7 



formation of low-metallicity stars to be at z>10 and its peak co-moving star-formation rate to 
be lower than 0.5 M Mpc~ 3 yr _1 . This upper limit is already of the same order of the peak 
star-formation rate of 0.2-0.6 M Q Mpc~ 3 yr _1 proposed by (47) and suggests that the peak 
star-formation rate might be much lower as proposed by (48). 
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Figure 1 Measurement, at the 68 % and 95 % confidence levels (including systematic uncertain- 
ties added in quadrature), of the opacity r 77 from the best fits to the Fermi data compared to 
predictions of EBL models. The plot shows the measurement at z^l which is the average red- 
shift of the most constraining redshift interval (i.e. 0.5<z<1.6). The Fermi-LAT measurement 
was derived combining the limits on the best-fit EBL models. The downward arrow represents 
the 95 % upper limit on the opacity at z=1.05 derived in (13). For clarity this figure shows only 
a selection of the models we tested while the full list is reported in Table SI. The EBL models 
of (49), which are not defined for E>250/ (1 + z) GeV and thus could not be used, are reported 
here for completeness. 
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Figure 2 Absorption feature present in the spectra of BL Lacertae objects as a function of 
increasing redshift (data points, from top to bottom). The dashed curves show the attenuation 
expected for the sample of sources by averaging, in each redshift and energy bin, the opacities 
of the sample (the model of (7) was used) and multiplying this average by the best-fit scaling 
parameter b obtained independently in each redshift interval. The vertical line shows the critical 
energy E crit below which <5 % of the source photons are absorbed by the EBL. The thin solid 
curve represents the best-fit model assuming that all the sources have an intrinsic exponential 
cut-off and that blazars follow the blazar sequence model of (32, 33). 
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Supplements 



Source Selection 

The second LAT AGN Catalog (16) contains many flat-spectrum radio quasars (FSRQs) and BL 
Lacertae (BL Lac) objects whose spectra extend significantly above lOGeV. The classification 
relies on the conventional definition of BL Lac objects outlined in (50-52) in which the equiv- 
alent width of the strongest optical emission line is <5 A and the optical spectrum shows a Ca 
II H/K break ratio C<0.4. Sources are then classified also according to the position of the peak 
of the synchrotron component as low- synchrotron-peaked (LSP, u peak <10 u Hz), intermediate- 
synchrotron-peaked (ISP, 10 14 <z/ peafc <10 15 Hz), and high- synchrotron-peaked (HSP, z/ peafe >10 15 Hz). 

FSRQs generally have soft spectra (photon index greater than 2.3) at GeV energies, making 
it more difficult to detect absorption features at energies greater than 10 GeV. Additionally 
their spectra might suffer from intrinsic absorption of gamma rays by the photons of the broad 
line region or of the accretion disk (17). This makes them non-ideal candidates to constrain 
the EBL. We therefore focused on the BL Lac objects. For computational reasons, only 50 
sources can be analyzed in one combined likelihood fit. In order to perform our analysis in 
three redshift intervals we selected the 150 BL Lacs that show the largest detection significance 
in the 3-10 GeV energy band (the significance is always larger than 3.5 a) with no selection on 
spectral shape. Most of the sources (100 out of 150) were detected with a significance >3a 
above 10 GeV already after two years of observations (21). In the second LAT AGN Catalog 
catalog (16) only ~190 BL Lac have a redshift measurement, so our sample is a representative 
set of Fermz-detected BL Lacs. Figure [ST] shows the redshift distribution of the BL Lac objects 
employed in this analysis. 

Data Analysis 

Each source in our sample is analyzed using 46 months^] of Fermi observations using version 
v9r27 of the Science Toolset The data were filtered, removing time periods in which the in- 
strument was not in sky-survey mode, and removing photons whose zenith angle is larger than 
100° to limit contamination from the Earth limb emission. We consider only photons collected 
within 15° of the source position with 1<E< 500 GeV. We employ the P7SOURCE_vdl in- 
strumental response function (IRF) and perform a binned likelihood analysis. The Galactic 
and isotropic diffuse emissions are modeled using respectively the gal_2yearp7v6_v0.fits and 
iso_p7v6source.txt templates. 

We rely on the Composite Likelihood tool of the Fermi software to perform the likelihood 
maximization. This allows us to fit simultaneously the data from different ROIs with the aim 
of constraining the scaling parameter b that is applied to the opacity curves predicted by each 

'From August 4th 2008 to June 1st 2012 
2 http://fermi. gsfc.nasa.gov/ssc/data/analysis/software/ 
3 Our analysis is robust against change of the dataset and IRF. 
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Figure SI Redshift distribution of the BL Lac objects used in this analysis. 

one of the EBL models we tested. The total number of free parameters across all ROIs in one 
redshift bin is approximately 1000. It is computationally unfeasible to fit all of these parameters 
simultaneously. We therefore proceeded in three steps. First we fit each ROI individually. 
The fit is performed on the entire energy band (1-500 GeV). The parameters (i.e. flux and 
photon index) of all the sources within 4° of the target source, along with the parameters of the 
diffuse components, are left free to vary. More distant sources have parameters frozen at the 
values measured in the second Fermi LAT catalog (27), unless the inspection of the residual 
map showed that a given source underwent strong variability. In that case, the normalizations 
of those sources were left free in the fit. The spectra of the sources of interest are modeled 
using a LogParabola model (21). We then proceed to re-optimize the parameters of the source 
of interest in each ROI up to that energy for which the EBL absorption becomes no longer 
negligible (see the definition of E crit in the main text). In the third step we fix the curvature of 
the LogParabola model to the value obtained in the previous fit. Finally, we use the Composite 
Likelihood to constrain the gamma-ray opacity. The Composite Likelihood allows the user to tie 
any parameters between any ROIs. In our case the only tied parameter is the renormalization 



15 



factor b of the opacity of a given EBL model. 

The significance of our finding can be evaluated using the Test Statistic (TS) as TS = 
2(logJtf(b) — logJ?(b = 0)), where J?f is the likelihood function and the b = case (i.e. no 
EBL absorption) represents the null hypothesis. Since the null hypothesis is a special case of 
the hypothesis we test, we expect the TS to be distributed as a x 2 with one degree of freedom 
(see also next sections). This allows one to transform the TS into the corresponding number of 
standard deviations of a Gaussian distribution as n a = y/TS. 

Modeling the Intrinsic Blazar Spectrum 

One of the basic assumptions of this analysis is that on average the spectrum of a BL Lac blazar 
can be adequately modeled using a LogParabola in the l-500GeV band. Figure [S2] shows 
the residuals of the best fits to all the sources in the z<0.2 interval. It is apparent that the 
LogParabola provides a good representation of the spectra of blazars in our sample. 

In an additional test we artificially decreased the critical energy E crit for all the sources in 
the z<0.2 interval from the typical > 120 GeV (for z<0.2) to ~40GeV which is representative 
of the z>0.5 case. Even in this case the results are unchanged, showing that the properties of 
the intrinsic spectrum can be determined over a more restricted energy range and can be safely 
extrapolated above the critical energy. 

Blazars are however known to have complex spectra that cannot always be modeled using 
simple empirical functions. To further test our basic assumption, we fit the spectra of bright 
blazars with coverage in the GeV and TeV band. We rely on the published spectra of Mrk 
421 (53), Mrk 501 (54) and RBS 0413 (55). We find that in the energy range relevant for the 
present analysis the LogParabola model provides an adequate fit to all these spectra. A further 
confirmation comes from our validation studies using simulations of broad-band blazar spectra 
(see next sections). We thus conclude that in a small energy range like the one adopted here 
(1-500 GeV), the LogParabola is an adequate representation of the intrinsic blazar spectrum. 

Results 

Table ISTl reports the results of the joint-likelihood fit for all the predictions of the EBL models 
we tested. In many cases, a given EBL model (because of tunable parameters or uncertainties 
connected to e.g, galaxy evolution, star-formation rate etc.) comprises several different predic- 
tions for the opacity. In this case, we use these prediction at face value and test them against our 
observations. All the models with the exclusion of rely on a standard concordance cos- 
mology (H =70km s _1 Mpc -1 , f2M=l-^A=0.3). For each model we report the significance of 
the absorption feature detected once we allow the opacity to be renormalized by b. The signifi- 
cances have been obtained as a square root of the sum of the TS in the 3 redshift intervals while 

4 The model of (22) adopts a value of Ho=65 km s _1 Mpc -1 . We checked that neglecting the change of Ho for 
the model of (22) introduces an error on the opacity of <7 % which has negligible influence on our analysis.. 
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Figure S2 Average of the residuals with respect to the best fit models for all the blazars in the 
z<0.2 redshift bin. 

the b parameter as reported in Table ISTIrepresents the weighted averagellof the 3 b values deter- 
mined in the 3 redshift bins. For example, for the model of (7) we measure: 6=1. 181^81 

(TS^4 

for z<0.2), 6=0.821^ (TS^7 for 0.2<z<0.5), and b=1.29t° H (TS^25 for 0.5<z<1.6). The 
weighted average yields (as reported in Table ISTI) b = 1.02 ± 0.23 and TS~36. 

We also quantify the level of compatibility of the prediction of a given model with the 
Fermi data, by comparing the likelihood of the original model (b = 1) to the one of the best 
fit scenario (6 free to vary). For b values significantly different than 1 a given EBL model 
predicts an attenuation larger than observed. Table [Si] shows that the intensity of the optical- 
UV background in 4 of the 15 models we tested is not compatible with the Fermi data at > 3 a. 

Figure [S3] shows that the significance of the absorption feature is detected collectively^ 
among the BL Lacs considered and is not attributable to just a few sources. Indeed, it is apparent 
that the TS increases linearly with the number of sources as expected in a background-limited 
scenario. 

As shown already in Figure [T] our analysis probes absorption signatures that are a factor 
~4 times weaker than those probed in (13) which relied on only one year of data. Even a 
conservative 95% upper limit derived from this analysis (b UL = 1.5 referred to the model 
of (7) <0.6 at z^l and E=77 GeV) is a factor ~2 below the one reported by (13). 

Our analysis relies on the assumption that BL Lac spectra (and HSP in particular) do not 
change dramatically across the z=0.2 redshift barrier. In the next section we address possible 
sources of systematic effects. 

5 The weights are where a\, is the uncertainty on the b parameters. 

6 In Figure[S3]the sources have been sorted in redshift from low to high redshift. 
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Table S 1 . Joint-likelihood results for different EBL models. 



MoHp1 a 


Ref b 




b d 


^li cttii fir* a Tir*^ of* n — 1 Rpipr^tion 


Stecker et al. (2006) - fast evolution 


(23) 


4.6 


0.10±0.02 


17.1 


Stecker et al. (2006) - baseline 


(23) 


4.6 


0.12±0.03 


15.1 


Kneiske et al. (2004) - high UV 


(22) 


5.1 


0.37±0.08 


5.9 


Kneiske et al. ( 2004 ) - best fit 


(22) 


5.8 


0.53±0.12 


3.2 


Gilmore et al. (2012) -fiducial 


(27) 


5.6 


0.67±0.14 


1.9 


Primack et al. (2005) 


(56) 


5.5 


0.77±0.15 


1.2 


Dominguez et al. ( 201 1 ) 


(25) 


5.9 


1.02±0.23 


1.1 


Finke et al. (2010) - model C 


(24) 


5.8 


0.86±0.23 


1.0 


Franceschini et al. (2008) 


(7) 


5.9 


1.02±0.23 


0.9 


Gilmore et al. (2012) - fixed 


(27) 


5.8 


1.02±0.22 


0.7 


Kneiske & Dole (2010) 


(26) 


5.7 


0.90±0.19 


0.6 


Gilmore et al. (2009) - fiducial 


(2) 


5.8 


0.99±0.22 


0.6 



a Models tested are implemented in the Fermi Science Tools. As an example the recent model of (49) which is not defined 
for E>250/(1 + z) GeV could not be used, but its predictions are for completeness reported in Figured! 

b Reference number in the 'References and Notes' section of the main text. 

c Significance, in units of cr, of the attenuation in the spectra of blazars when a given EBL model is scaled by the factor b. In 
this case 6=0 (i.e. no EBL absorption) constitutes the null hypothesis. 

d This column lists the maximum likelihood values and 1 a confidence ranges for the opacity scaling factor. 

c Here the b=l case (i.e. EBL absorption as predicted by a given EBL model) constitutes the null hypothesis. This column 
shows the compatibility (expressed in units of a) of the predictions of EBL models with the Fermi observations. Large values 
mean less likely to be compatible. 




Number of Sources 



Figure S3 Increase in the TS value of the (renormalized) EBL model of (7) produced in the 
joint-likelihood fit (to the 0.5<z<1.6 interval) while adding one source at a time. The sources 
have been sorted in redshift (from lowest to highest). The dashed line shows the best-fit lin- 
ear increase of the TS with the number of sources. The inset shows the best-fit value of the 
renormalization parameter b applied to the opacity predicted by (7) (see text for details). 

Validation of the Analysis using Simulations 

We validated our analysis using simulated datasets starting from physically motivated spectral 
energy distributions (SEDs) of BL Lacs. We rely on SEDs that reproduce well the range of 
peak curvatures, peak frequencies (for both the synchrotron and high-energy component), and 
gamma-ray photon indices observed for the LAT-detected BL Lacs. These SEDs have been 
produced using the same numerical code presented in (57) and take into account all the im- 
portant effects that contribute to determine the intrinsic curvature of the gamma-ray spectrum 
of the Fermi blazars. Of particular importance are those effects that depend on the curvature 
of the electron distribution, as well as those related to the Thomson (TH) to Klein-Nishina 
(KN) transition, in the inverse Compton process (IC). The latter effect is crucial since it is well 
known (57, 58) that the TH/KN transition in the IC cross section can result in a steepening of 
the high-energy branch of the IC SED, compared to the low-energy one. Our simulations fully 
take into account this effect, allowing us to study the potential bias in the EBL estimate. 

We performed ~500 simulations randomly selecting 50 SEDs from our template library. 
To each SED we assign a redshift (in the 0.5<z<1.6 range) and a 0.1-100 GeV flux randomly 
extracted from those of BL Lacs in the second LAT AGN Catalog (16). We then proceeded to 
simulate the data expected from those sources based on Fermi's instrument response function 
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and pointing history. Each one of those realizations was processed with the same analysis chain 
used on actual data. 

Figure [S4] presents the distribution of TS (left panel) and best-fit b parameters (right panel) 
for the simulated blazar population under the following two scenarios. First, no EBL attenuation 
was included, in order to check whether our analysis would pick a signal consistent with EBL 
attenuation when this was not present (i.e. a false positive). The dashed lines in the panels show 
that our analysis yields very small TS values and b values compatible with zero, as expected for 
a robust analysis. 

In the second scenario, the SED models were attenuated using the opacities of the model 
of (7). The results summarized in Figure IS4l show that TS values >20 and b values compatible 
with 1 (within the statistical uncertainty) are correctly retrieved. Moreover, the simulations 
show that our analysis is free of any major systematic uncertainty since the average values of b 
retrieved for both scenarios are compatible with the expected ones (0 and 1 respectively) within 
1%. 




Test Statistics scaling parameter b 



Figure S4 Results from joint-likelihood fits to the simulated datasets. The left panel shows the 
TS of the detection of the attenuation produced by the EBL for the case in which no EBL was 
applied to the SEDs (dashed line) and the case for which an attenuation consistent with the one 
of (7) was used (solid line). The right panel shows the best-fit b parameter for the same two 
scenarios. 
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Additional Tests 



A number of additional tests have been performed to test the stability of the results presented in 
the previous sections. These are documented below. 

Sample Selection 

The selection criteria used in the previous sections were chosen a priori to yield a uniform 
source population that could be used to probe in a sensitive and least-biased way the EBL- 
induced absorption feature. Several tests have been performed to check this assumption. 

The BL Lac population is known to evolve in redshift with HSP objects detected predom- 
inantly at low redshift and LSPs detected predominantly at larger redshifts with ISPs bridging 
the gafl This could artificially induce a decline in the combined (produced by all subtypes) 
BL Lac spectrum at an energy that decreases with increasing redshift, and hence could mimic 
the effect of EBL-caused absorption. To test whether such bias affects our results we run the 
analysis separately for the three source classes (HSPs, ISPs and LSPs). For each of them we 
define two redshift bins: z<0.2 and z>0.2. In the lowest redshift bin, due to the low statistics^ 
the joint likelihood yields no significant detection of EBL softening. For the HSPs in the low 
redshift bin, the joint likelihood yields 6=1. 46^' and a TS=4.1. At higher redshift, the results 
are: b =1.04+°;!| (TS=24), b =2.13^ (TS=9), and b =0.43lg;?J (TS=1), for HSPs, ISPs and 
LSPs respectively. The weighted average of the above values yields b =1.04+° 24 and a total 
TS of ~38. These results are in very good agreement with the ones presented in the previous 
sections and show that high-redshift LSPs are not responsible for the observed spectral soften- 
ing, and HSPs are mainly responsible for it. Moreover, since the signal is dominated by only 
one source population (namely HSPs) the probability that selection effects from changing the 
sample composition with redshift strongly affect our result is rather low. 

We also checked if the highest redshift HSPs are driving the result. To this extent we isolated 
(from the above sample) the 6 HSPs with the highest redshift. The joint likelihood yields a TS=6 
and b =1.991^^ indicating that the detection of the EBL cut-off is not caused by a few high 
redshift sources. 

We also tested the population of BL Lacs that did not pass our selection criteria. This 
comprises 32 objects detected with a rather low significance in the 3-10 GeV over a redshift 
range 0.03-1.9 (average of 0.4). The joint likelihood yields b=l.93t.fH and TS=3.5 compatible 
with our result. Thus, including low significance sources marginally improves the results of our 
analysis. 

7 The sample used in the previous sections contains for z<0.2 35 HSPs, 10 ISPs and 5 LSPs. The 0.2<z<0.5 
sample comprises 27 HSPs, 18 ISPs, 5 LSPs, while the z>0.5 sample contains 10 HSPs, 19 ISPs, and 21 LSPs. 

8 In the z<0.2 bin there are respectively 41 HSPs, 10 ISPs and 5 LSPs; while in the z>0.2 there are 50 HSPs, 
30 ISPs, and 23 LSPs. 
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Binning in Redshift 

We bin the sources in redshift to cope with the maximum number of free parameters (i.e. 100) 
that can be optimized in a single maximum likelihood fit. Our analysis is however independent 
of the redshift binning used. This is already evident from the tests reported above using sep- 
arately HSPs, ISPs and LSPs and a different redshift binning (with respect to the one used in 
the main analysis). However, we performed an additional test and divided the most constrain- 
ing redshift bin (i.e. 0.5<z<1.6) into halves (at z^0.75) containing 25 sources each. When 
fitting the Fermi data using the model of (7), we obtain 6=1.7 11^82 (TS=10) and 6=1.17+0 39 
(TS=15.4) for the 0.5<z<0.75 and 0.75<z<1.6 bins respectively. The weighted average of the 
above values yields 6=1.26+0 35 an( ^ me total TS=25.4. This can be compared to 6=1.29+q 3q 
and TS=25.1 obtained for the full 0.5<z<1.6 redshift bin. Our analysis is thus robust against 
choices of the redshift bins. 

Source Variability 

Blazars are inherently variable objects with variability in flux of up to a factor 10 or more. 
Throughout this work we use spectra of blazars accumulated over 46 months of observations 
and blazar variability might in principle bias our result. In the Fermi sample, the most variable 
BL Lac objects are those belonging to the LSP and ISP class and it has been shown already that 
the entire analysis can be confirmed without using those two classes. 

To test possible biases deriving from variability we selected the 30 most variable, according 
to (16, 21), and most significant (in the 3-10 GeV band) BL Lac objects. Repeating the entire 
analysis using this sample, we obtain b = I.2O+0 52 and a TS=9. These results are in good 
agreement with the ones presented in the previous section and show that blazar variability does 
not affect our results. 

Critical Energy 

Our measurement relies on the extrapolation of the intrinsic (i.e. unabsorbed) spectrum above 
the critical energy E crit . Aggressive choices of this energy (i.e. large values) might bias the 
result towards a low level of opacity. We show that this is not the case for our measurement and 
that our analysis is robust against conservative choices of E crit . In order to demonstrate it, we 
imposed E crit =l0 GeV for any source in the most constraining redshift interval: i.e. 0.5<z< 1 .6. 
For most EBL models, with the exception of those presented in (23), the opacity is negligible 
at 10 GeV up to z^l.6. When fitting the Fermi data using the model of (7), our analysis yields 
a value of &=1. 24+^37 and TS?^22 for the constant E crit =\0GeV case. This can be compared 
to 6=1.29+Q 3g and TS?^25.1 for the variable E crit case presented in the main text. We thus 
conclude that our analysis is robust against any conservative choice of E crit . 
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Systematic Uncertainties of the Instrument Response Function 

In order to gauge the systematic uncertainties of this analysis we use the IRF bracketing method 
as described in (59). By deriving two different sets of IRFs and repeating the entire analysis we 
find that the systematic uncertainty on the opacity r 77 is of the order ~7 %. 

Moreover, the entire analysis has been repeated using the P7CLEAN_V6 photon selection 
and IRF. The results presented in the above sections are fully confirmed and the systematic 
uncertainty on the opacity r 77 due to changing event selection and IRF is of the order 3 %. 
Thus, we consider the total systematic uncertainty on r 77 to be ~10 %. 

Source Intrinsic Effects 

The blazar evolution with luminosity as described by the blazar sequence (31-34) and the 
change of blazar type with redshift might in principle produce an intrinsic cut-off that changes 
with redshift. An intrinsic spectral cut-off might be detected if the part of the blazar spectrum 
sampled by Fermi corresponds to the tail of the electron distribution. In our case, this might 
happen only for the LSP sources and it has already been shown (see previous section) that their 
contribution to the total signal is negligible. On the other hand most of the signal is dominated 
by HSPs that are detected by Fermi right below or at the peak of the inverse Compton com- 
ponent, thus excluding the possibility that their emission is produced in the tail of the electron 
distribution. 

However, in order to study the compatibility of the detected signal with an intrinsic origin, 
we assumed that all blazars have an exponential cut-off at a source-frame energy E . We also 
assumed that the blazar spectrum moves to lower frequency for increasing source luminosity 
as dictated by the so-called blazar sequence (31-34). This might represent the case in which 
the maximum electron energy depends on the jet power or the luminosity. Figure [2] shows that 
this model (E is fitted to the data) provides a poor representation of the Fermi data. Even 
assuming that the maximum electron energy depends on a variable power of the luminosity 
does not improve the results. The result is even worse if the blazar sequence is not invoked. We 
thus conclude that the impact of source intrinsic effects on our analysis is likely minor. 

Reprocessed Emission 

The electron-positron pairs generated in the interaction of gamma rays and EBL photons can 
initiate a cascade by subsequent Compton scattering of the photons of the Cosmic Microwave 
Background. Typically the cascades originated by TeV gamma rays are reprocessed in the 
GeV energy range. In the case of a weak intergalactic magnetic field (IGMF), the pairs are 
not deflected out of the line of sight (60) and the reprocessed emission can be a substantial 
fraction of (or even dominate) the total source signal in the <100 GeV band. The intensity of the 
reprocessed emission depends primarily on the EBL density, on the intensity of the IGMF, on its 
coherent length and on the position of the peak of the high-energy component of a blazar (61). 
For typical coherent lengths of ~ 1 Mpc and SED peaks located at or below 10 TeV, lower limits 
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of B>10~ 15 G have been inferred (61, 62). For these intensities of the IGMF the reprocessed 
emission is expected to be subdominant with respect to the intrinsic component (63). While 
it is true that these estimates are based on an implied activity of the source of a few million 
years (64), the majority of the Fermi BL Lacs is expected to have the peak of the high-energy 
component located at < 1 TeV (65) greatly reducing the amount of reprocessed emission in the 
GeV band (61). Similar considerations (i.e. suppression of the reprocessed emission by the 
IGMF) hold also for the case in which blazars are also sources of ultrahigh energy cosmic 
rays (66-68). 
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